#include<stdio.h>
#include<stdlib.h>
#include<math.h>

/* Função para calcular o tempo de execução */
double tempo()
{
    struct timeval tv;
    gettimeofday(&tv,0);
    return tv.tv_sec + tv.tv_usec/1e6;
}

/* Inicializa uma matriz */
double **initMatriz(int ordem){
  int i;
  
  /* Cria e aloca espaço para a matriz */
  double **matrix;
  matrix = malloc(ordem*sizeof(double *));
  for (i = 0; i < ordem; i++)
    matrix[i] = malloc(ordem*sizeof(double)); 
  
  return matrix;
}

/* Inicializa um vetor */
double *initVector(int ordem){
  /* Cria e aloca espaço para o vetor */
  double *vector;
  vector = malloc(ordem*sizeof(double));
  
  return vector;
}

/* Libera os espaços em memoria alocados pelas matrizes */
void free_matrix(double **mat, int rows){
    int i;
    for(i=0; i<rows; i++){
        free(mat[i]);
    }
    free(mat);
}

/* Multiplica uma Matriz por um Vetor */
void multiMat_Vec(double *result, double **mat, double *vec, int ordem){
  int i, j;
  double soma;
  
  for(i = 0; i < ordem; i++){
    soma = 0;
    for(j = 0; j < ordem; j++){
      soma += mat[i][j] * vec[j];
    } 
    result[i] = soma;
  }  
}

/* Verifica o erro */
int verificaErro(double *X, double **MA, double *MB, int ordem, int linha_teste, double erro){
  int i;
  double result = 0;
  
  for(i = 0; i < ordem; i++){
    result += MA[linha_teste][i] * X[i];
  }
  
  if(fabs(result - MB[linha_teste]) < erro) return 0;
    
  return 1;
}




int main(int argc, char *argv[]){
  int i, j, iteracao;
  double result;
  
  /* Variáveis para registrar o tempo de início e término de execução */
  double start, end;
  
  /* Ordem da matriz */
  int J_ORDER;
  
  /* Fila a ser avaliada */
  int J_ROW_TEST;
  
  /* Erro permitido */
  double J_ERROR;
  
  /* Número máximo de iterações */
  int J_ITE_MAX;
  
  /* Matrizes utilizadas */
  double **MA, *MB, **LR, *b;
  
  /* Vetor de resposta */
  double *X, *Xaux;
  
  /* Ponteiro para o arquivo de entrada e arquivos de saída (estatísticas e vetor X resposta) */
  FILE *fp;
  
  /* Registra tempo de início */
  start = tempo();
  
  /* Abre o arquivo, de acordo com o parâmetro de entrada */
  if(argc <= 1){
    printf("Numero insuficiente de parametros. Informe o nome do arquivo.\n");
    return -1;
  }
  fp = fopen(argv[1], "r");
  if(fp < 0){
    printf("Arquivo não encontrado.\n");
    return -1;
  }
  
  /* Lê os 4 primeiros parâmetros */
  fscanf(fp, "%d%d%lf%d", &J_ORDER, &J_ROW_TEST, &J_ERROR, &J_ITE_MAX);
  
  /* Cria a matriz A, com o parâmetro fornecido */
  MA = initMatriz(J_ORDER);
  
  /* Lê matriz MA (ainda não otimizado) */
  for( i = 0; i < J_ORDER; i++ ){
    for( j = 0; j < J_ORDER; j++ ){
      fscanf(fp, "%lf", &(MA[i][j]) );
    }
  }
  
  /* Cria o vetor B, com o parâmetro fornecido */
  MB = initVector(J_ORDER);
  
  /* Lê vetor MB (ainda não otimizado) */
  for( i = 0; i < J_ORDER; i++ ){
      fscanf(fp, "%lf", &(MB[i]) );
  }
  
  /* Fecha o arquivo */
  fclose(fp);
  
  /* Passos para a matriz LR, que já é a soma da matriz L com a matriz R, invertendo o sinal (simplificação do algoritmo explicada no relatório) */
  LR = initMatriz(J_ORDER);
  for( i = 0; i < J_ORDER; i++ ){
    for( j = 0; j < J_ORDER; j++ ){
      if(i > j) LR[i][j] =  - MA[i][j]/MA[i][i]; /* Matriz L */
      else if(i < j) LR[i][j] = - MA[i][j]/MA[i][i]; /* Matriz R */
      else LR[i][j] = 0.0; /* Diagonal zero */
    }
  }
  
  /* Cria o vetor b, que é são os elementos de MD divididos pelos elementos da diagonal de MA */
  b = initVector(J_ORDER);
  for( i = 0; i < J_ORDER; i++ ){
      b[i] = MB[i]/MA[i][i];
  }
  
  /* Prepara o vetor x */
  X = initVector(J_ORDER);
  /* Inicializa com valores iniciais. No caso, coloquei zeros */
  for( i = 0; i < J_ORDER; i++ ){
    X[i] = 0.0;
  }
  
  /* Vetor auxiliar para computar o resultado intermediário de LR*X */
  Xaux = initVector(J_ORDER);
  
  
  /* Começa a iteração do método */
  for(iteracao = 0; iteracao < J_ITE_MAX; iteracao++){
    
    /* Faz a operação -(L + R)*X{k} */
    multiMat_Vec(Xaux, LR, X, J_ORDER);
    
    /* Calcula o X{k+1} */
    for(i = 0; i < J_ORDER; i++){
      X[i] = Xaux[i] + b[i];
    }
    
    /* Se está dentro do erro permitido, termina a execução */
    if (verificaErro(X, MA, MB, J_ORDER, J_ROW_TEST, J_ERROR) == 0){
      break;
    }
  }
  
  /* Registra tempo de fim */
  end = tempo();
    
  /* Imprime os resultados no terminal */
  printf("----------------------------------------\n");
  printf("Iterations = %d\n", iteracao);
  
  result = 0;
  for(i = 0; i < J_ORDER; i++){
    result += MA[J_ROW_TEST][i] * X[i];
  }
  printf("RowTest: 2 => [%lf] =? %lf\n", result, MB[J_ROW_TEST]);  
  printf("Time = %f seconds\n", end-start);
  printf("----------------------------------------\n");
  
  /*=(Adicional para ver os resultados)=======================================*/
  /* Registra a execução em um arquivo de estatísticas */
  fp = fopen("statistics.txt", "a");
  if(fp < 0){
    printf("Nao foi possivel adicionar resultado ao arquivo de estatisticas.\n");
    return -1;
  }
  fprintf(fp, "----------------------------------------\n");
  fprintf(fp, "Method: Sequential\n");
  fprintf(fp, "Matriz order: %d\n", J_ORDER);
  fprintf(fp, "Iterations = %d\n", iteracao);
  fprintf(fp, "RowTest: 2 => [%lf] =? %lf\n", result, MB[J_ROW_TEST]);  
  fprintf(fp, "Time = %f seconds\n", end-start);
  fclose(fp);
  
  /* Agora registra o vetor X em um arquivo específico (este não será incremental) */
  fp = fopen("resultX.txt", "w");
  if(fp < 0){
    printf("Nao foi possivel criar um arquivo para armazenar o vetor X.\n");
    return -1;
  }
  fprintf(fp, "vectorX = {\n");
  for( i = 0; i < J_ORDER; i++ ){
    fprintf(fp, "      %lf,\n", X[i]);
  }
  fprintf(fp, "}\n");
  fclose(fp);
  /*======================================================================*/
  
  /* Libera da memória as matrizes e vetores usados */
  free_matrix(MA, J_ORDER);
  free_matrix(LR, J_ORDER);
  free(MB);
  free(b);
  free(X);
  free(Xaux);
    
  return 0;
}
